import numpy as np
from scipy import stats

def find_min_n_for_acceptance(p0=0.1, alpha=0.1, delta_start=0.001, delta_end=0.05, delta_step=0.001, max_n=100000):
    """
    在90%置信度下，找到满足接收条件的最小样本量
    接收条件: p̂ + z_{α/2} · SE ≤ p₀
    """
    z_alpha_2 = stats.norm.ppf(1 - alpha/2)
    print(f"z_(α/2) = {z_alpha_2:.4f}")
    print(f"置信水平 = {90}%")
    print(f"标称次品率 p₀ = {p0}")
    print(f"接收条件: p̂ + {z_alpha_2:.4f} · SE ≤ {p0}")
    print(f"{'p':>8} {'min_n':>8} {'示例k':>8} {'p̂':>8} {'SE':>10} {'接收条件':>12} {'满足条件':>8}")

    delta = delta_start
    while delta <= delta_end:
        p = p0 - delta  # 实际次品率低于标称值
        n = 10
        found = False
        while n <= max_n:
            # # k ~ B(n, p)
            # k = np.random.binomial(n, p)
            # # 样本次品率
            # p_hat = k / n
            k = n * p
            p_hat = p
            # 标准误差 - 处理边界情况
            if p_hat == 0:
                # 当p̂=0时，使用保守估计
                SE = np.sqrt(0.5 * 0.5 / n)  # 使用0.5作为保守估计
            elif p_hat == 1:
                # 当p̂=1时，使用保守估计
                SE = np.sqrt(0.5 * 0.5 / n)  # 使用0.5作为保守估计
            else:
                SE = np.sqrt(p_hat * (1 - p_hat) / n)
            
            # 接收条件: p̂ + z_{α/2} · SE ≤ p₀
            acceptance_criterion = p_hat + z_alpha_2 * SE
            # 检查是否满足接收条件
            if acceptance_criterion <= p0:
                print(f"{p:8.4f} {n:8d} {k:8.0f} {p_hat:8.4f} {SE:10.6f} {acceptance_criterion:12.6f} {'是':>8}")
                found = True
                break
            n += 1
        if not found:
            print(f"{p:8.4f} {'未找到':>8}")
        delta += delta_step

if __name__ == "__main__":
    print("=== 90%置信度接收条件分析 ===")
    # 90%置信度下的接收条件
    find_min_n_for_acceptance(p0=0.1, alpha=0.1, delta_start=0.001, delta_end=0.051, delta_step=0.0001, max_n=10000)
